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ABSTRACT 

The interaction of ionizing and far-ultraviolet radiation with the interstellar medium is of great 
importance. It results in the formation of regions in which the gas is ionized, beyond which 
are photodissociation regions (PDRs) in which the gas transitions to its atomic and molecular 
form. Several numerical codes have been implemented to study these two main phases of the 
interstellar medium either dynamically or chemically. In this paper we present TORUS-3DPDR, 
a new self-consistent code for treating the chemistry of three-dimensional photoionization and 
photodissociation regions. It is an integrated code coupling the two codes TORUS, a hydrody¬ 
namics and Monte Carlo radiation transport code, and 3D-PDR, a photodissociation regions 
code. The new code uses a Monte Carlo radiative transfer scheme to account for the prop¬ 
agation of the ionizing radiation including the diffusive component as well as a ray-tracing 
scheme based on the HEALPIX package in order to account for the escape probability and 
column density calculations. Here, we present the numerical techniques we followed and we 
show the capabilities of the new code in modelling three-dimensional objects including sin¬ 
gle or multiple sources. We discuss the effects introduced by the diffusive component of the 
UV field in determining the thermal balance of PDRs as well as the effects introduced by a 
multiple sources treatment of the radiation field. We find that diffuse radiation can positively 
contribute to the formation of CO. With this new code, three-dimensional synthetic obser¬ 
vations for the major cooling lines are possible, for making feasible a detailed comparison 
between hydrodynamical simulations and observations. 
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1 INTRODUCTION 

The interstellar medium (ISM) consists predominantly of hydrogen 
(70%), followed by helium (28%) and heavier elements (2%). Its 
interaction with the ultraviolet (UV) radiation emitted by massive 
stars is of great interest. Extreme UV photons carrying more en¬ 
ergy than the ionization potential of hydrogen (13.6 eV) ionize the 
gas. Any excess over the ionization potential is then transformed 
into kinetic energy of the liberated electrons (photoelectrons). Re¬ 
gions containing a signifi cantly high fracti on of photoelectrons are 
known as ‘H II regions’ ( IStr6merenlll939tl and have temperatures 
of the order of 5, 000 — 15, 000 K which is higher than t hat of 
the surrounding gas (^ 100 K) resulting in their expansion dKahnl 
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Il954l ; l^itzeilll978h . The edge of the ionized medium is known as 
an ‘ionization front’ and separates the H II region from the rest of 
the ISM. The ionization front can be considered as a sharp discon¬ 
tinuity in which the degree of ionization and the gas temperature 
drops abruptly, over a few UV photon mean free paths. The ma¬ 
terial ahead of this front is dominated by far-UV radiation (FUV) 
carrying energy Q < hv < 13.6 eV. Photons in this energy regime 
can photodissociate CO and H 2 molecules creating so called ‘pho¬ 
todissociation region£_(PDRs]_^nown_also as ‘photon-dominated 
regions’, see iHollenbach & Tielenslll999L for a review). Further 
ahead the radiation field is attenuated to the extent that it does not 
affect quiescent molecular zones. 

Studies of the different phases of the ISM are of great inter¬ 
est as they can lead to an in-depth understanding of the mecha¬ 
nisms responsible for its structure and dynamical evolution. This in 
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turn may answer important questions concerning the contribution 
of s tellar feedback in driving turbulence at small (subpc) scales (see 
e.g. lKlessen & Gloveill20l4 . for a more general discussion on ISM 
turbulence) well as the modes that trigger the formation of new 
stars beyond the expanding ionized regions. Studying PDRs is key 
to understanding the role played by the FUV photons in governing 
the physical and chemical structure and determining the thermal 
balance of the neutral ISM in galaxies. 

Several groups worldwide have made significant efforts in 
simulating the (hydro-)dynamical evolution of the ISM using 
different computational techniques such as Smoothed Parti¬ 
cle Hydrodynamics (SPH) and Adaptive Mesh Refinement 
(AMR). These have led to the modelling of star formation 
processes, as well as feedback from massive stars and mag- 
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of data from the atomic and molecular lines emitted at different 
wavelengths by the different phases of the ISM. However, this does 
not directly reveal the exact three-dimensional density and velocity 
structure of the observed objects. There is thus the need for a tool 
which acts as a “common interface” between the simulated hydro- 
dynamical ISM structure and the observed line emission. Towards 
this goal, efforts have been made to model the complex chemical 
processes occurring in the different phases of the ISM of particular 
relevance to PDRs. Coupling radiation transport techniques with 
hydrodynamics and the gas chemistry, however, is complicated 
and poorly developed as yet, al though significant e f forts h ave 
been made in this direction by 


ly Glover & Mac Lo^ ll2007a|[9): 

IdTOlOt) andlLevrier et al .U2012I). 


iDobbs et al.l ( l2008l) : lGlover et ^ 

Numerical codes dedicated to modelling the chemistry 
in photoionized regions and PDRs have been presented in 
the past. Codes fo r photoionized regions include CLOUDY 
dperland et al.l Il998h while a three-dime nsional treatment has 
been achieved in the MO C ASSIN code l Ercolano_etalJ_ 200^ 
Wood. Mathis. & Ercoland |2()04 lErcolano. Barlow. & Storey 


2005) and in the TORUS code l HarriesI 200d : Haworth & Harries 


2012h which are based on Monte Carlo (MC) methods first 


described by ILucvI dl999l) . On the other hand, several codes 
treating PDRs have been pre sented in the past decade or so . 
These include UCL JDR dPanadopo ulos. Th i. & Vitil l2002l: 
Bell et al. 200^ 20061). C LOUDY dFerlandelj^ 199i: Abel et all 
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Shaw et akl 


COSTAR 
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iRbllig et"^ d2007h provide an interesting review of many of 
those PDR codes as well as benchmarking tests pinpointing the 
differences. The codes are able to treat complicated chemical net¬ 
works and thus obtain the abundances of many chemical species, 
the line emissivities of different coolants, as well as temperature 
profiles within a given PDR. While most attention has been paid to 
an accurate chemical PDR modelling, th e majority of PDR codes 
treat one-dimensional structures. Recentlv lBisbas et al.l ( l2012h have 


implemented the three-dimensional code 3D-PDR which is able 
to tr eat PDRs of any given arbitrary density d istribution. In addi¬ 
tion lAndree-Labsch. Ossenkopf. Rollid ( l20I4h also improved the 
KOSMA-r PDR code to handle 3D pixels (voxels) which mimic the 
fractal structure of the ISM. These new implementations meet the 
need for PDR codes to treat complicated geometrical structures, 
since the chemistry can depend on it. 

In this paper, we present the first self-consistent fully three- 
dimensional unified code for treating photoionized and photodisso- 
ciated regions simultaneously with arbitrary geometrical and den¬ 
sity distributions. This new code is able to calculate i) the struc¬ 
ture of circumstellar and interstellar radiation fields using an MC 
method that treats the diffuse component and ii) a detailed PDR 
chemistry using a HEALPIX ray-tracing scheme to estimate the 
cooling and heating rates as well as line emissivities, the abun¬ 
dances of chemical species and temperature profiles, following the 
methods of the 3D-PDR code. It therefore models the wide range 
of observable chemical processes which take place in the transition 
zones from ionized to molecular regions. 

Our paper is organized as follows. In Sections|2and[^we give 
an overview of the 3D-PDR and TORUS codes respectively. In Sec- 
tion|4]we show the method we followed for coupling the two codes. 
In Section [3 we benchmark our new code while in Section]^ we 
show the capabilities of the code in simulating three-dimensional 
structures. We conclude in Section|7]with the description of the ef¬ 
fects introduced in simulating photoionized and photodissociation 
regions in full three dimensions. 


2 THE 3D-PDR CODE 

The 3D-PDR code dBisbas et all |2012| . hereafter ‘B12’) is a 
three-dimensional time-dependent astrochemistry code which has 
been designed to treat PDRs of arbitrary density distribution. It 
is a further development of the one-dimens ional UCL _ pdr c ode 
and it uses the chemical model features of iBell et aO ( I 2 OO 6 I) . It 
solves the chemistry and the thermal balance self-consistently 
in each computational element of a given cloud. 3D-PDR has 
been fully b enchmarked agains t other PDR codes according to 
the tests of iRollig et alj (l2007h and has b een used in variou s 
applications e.g. lO ffner et alJ (1201.4 l2014l): iB isbas et all ll2014l) : 
iBisbas. Papadopoulos. & Vitil ( l2015h : lGaches et al.l(l2015h . 

We list below a brief overview emphasizing on the ray-tracing 
and UV treatment as well as the model chemistry used in 3D-PDR. 
For full details about the code see B12. 


2.1 Ray-tracing and UV treatment 

3D-PDR uses a r ay-tracing scheme based on the HEALPIX 
dGorski et alj|200^ package. Along each ray, the elements clos¬ 
est to the line-of-sight are projected creating a set of points called 
‘evaluation points’. The properties of the evaluation points are iden¬ 
tical to those of the projected elements. With this technique we are 
able to calculate: i) the column densities of species for a random 
element along a particular direction, ii) the attenuation of the far- 
ultraviolet radiation in the PDR, and iii) the propagation of the far- 
infrared/submm line emission out of the PDR. 

The treatment of the UV radiation in the 3D-PDR code is ini¬ 
tially estimated by invoking the on-the-spot approximation i.e. by 
neglecting the diffuse component of the radiation field. The attenu¬ 
ation of the UV field, x, at a randomly given point, p, is evaluated 
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using the equation 

= ( 1 ) 

^ q=l 

where Xo(q) is the magnitude of the unattenuated field strength 
along the HEALPIX direction q, Me = 12 x 4^ is the number of 
HEALPIX rays at level of refinement i, ruv = 3.02 is a dimen- 
sionless factor c onverting the visual extinction to UV attenuation 
( Ib^ ^t_aj|200^), and ./dv is the visual extinction along q defined 
as 

^v(q) = 4lv,o [ nudr . (2) 

Jo 

In the above equation Av,o = 5.3 x mag ctM and the inte¬ 

gration corresponds to the column of the H-nucleus number density 
riH along the ray of length L. We note that TORUS-3DPDR performs 
radiative transfer with frequency dependent opacity, however for 
the purposes of this initial comparison in the present work and con¬ 
sistency with the major ity of the PDR co des we have restricted ruv 
to the above value (c.f. IPinte et aLll2009h . 

Equation [T] has been extensively used in several codes treat¬ 
ing PDRs. As described in iRbllig et alJ l l2007h . most of the models 
adopt a plane-parallel geometry illuminated from one or from both 
sides. However, the majority of them do not account for the geo¬ 
metrical dilution i.e. that a spherically emitted radiation field de¬ 
creases in intensity with the square of the distance. Accounting for 
the latter, the attenuation along a HEALPIX direction q is expressed 
as follows: 

X(Aq)^Xo(q)e — -v(q)_^, (3) 

where is the position of the ionization front from the excit¬ 
ing source, and r is the distance of a particular point p from the 
ionization front and inside the PDR (see Appendix for the cor¬ 
responding derivation). As we explain below, this is an important 
factor especially when a PDR model is being used to reproduce the 
observed data of an object. In TORUS-3DPDR this attenuation natu¬ 
rally comes from the Monte Carlo (MC) radiative transfer treatment 
and this is what is being used throughout this paper unless other¬ 
wise stated. 


of molecular hydro gen formation on dust grains i s calculated us¬ 
ing the treatment of Cazaux & TielensI | |2002L[200^ while the ther¬ 
mally averaged stic king coefficient of hy d rogen atoms on dust 
grains is taken from iHollenbach & McKea ( Il979l) . The dust tem¬ 
perature at each poin t in the density distribution is calculated us¬ 
ing the treatment of IHollenbach. Takahashi. & Tielen^ l ll99lh to 
account for the grain heating due to the incident FUV photons. Fu- 
ture updates will include mo re de tailed heating functions such as 
IWeingartner & Draind ( 1200 ll) and lComniegne et"^ (1201 ll) for the 
grain charge photo-electric heating and grain temperature respec¬ 
tively. 


3 THE TORUS CODE 


TORUS is a grid-based three-dimensional radiation transport and 
Eulerian hydrodynamics code that uses an octree adaptive (non- 
uniform) grid. 

It was initially designed to perform three-dimensional cal¬ 
culations, including the treatment of polarization and Mie or 
Rayleigh scattering. In this form it was used to analyze syn¬ 
thetic spectral line observations of stellar winds that are ro- 
tationall v distorted by rapid stellar rotation or which contain 
clumps dHarriesluOOOl) . With the addition of a dust treatment it 
was also used to model observations of the Wolf-Rayet (WR) 
star binary WR137 in an effort to provide an explanation fo r 
observed polarization variability ( iHmies. Babler. & FoxI I2OOOI1 . 
It has since continued developing and been applied to mod¬ 
els of accretion on to T Tauri stars dvink, Harries. & Drew! 


I 2 OO 5 I: ISvmington. Harries. & Kurosawal 20051). discs around Her- 


big AeBe stars dTannirkulam et al .| 120081), Raman-scattere d line 
formation in symbiotic binaries ( Harries & Howard] 1 1997h . dust 
emission and molecular line formation i n star forming regions 
dKurosawa et al.ll200l: iRund le et al.ll 201(lll and synthetic galactic 
HI observations ( Acreman et al. 2010h . Most recently it has been 
used in radiatio n hydrodynamic applications, both by coupling it 
to an SPH code dAcreman. Harries. & Rundi3l2010ll and in a self- 
contained mann er by performing the hydro dynamics calculation on 
the TORUS grid dHaworth & Harrie3l2012h . 

We give below a brief overview emphasizing on the key fea¬ 
tures of the MC photoionization routines. 


2.2 Model chemistry 


3.1 Monte Carlo Photoionization 


The code uses the most recent UMIST 2012 chemical network 
database dMcElrov et al.l 1201 3h . This network consists of 215 
species and more than 3000 reactions. However, in this paper we 
will consider only a subset of this network consisting of 33 species 
(including e“) and 330 reactions. We make use of the SUNDIALS 
packag^in order to construct the appropriate set of ordinary differ¬ 
ential equations (DDEs) for describing the formation and destruc¬ 
tion of each species as well as the associated Jacobian matrices to 
compute the abundance of each chemical species. 

For the particular cases of H 2 a nd CO photodis socia¬ 
tion rates, we adopt the tr eatments of lUee et al] dl996ll and 
Ivan Dishoeck & BlackI dl988tl and we use the additional tabulated 
shielding functions provided in those pa pers. We also account for 
the shielding of C I using the treatment of lKamp & Bertoldil ( l2000tl 
in order to estimate the photoionization rate of carbon. The rate 


^ http://computation.llnl.gov/casc/sundials/main.html 


TORUS performs photoionization calculations using an iterative 
M C photon energy pack et pr opagating routine, similar t o tha t 
of lErcolano et al. ( l2003h and IWood. Mathis. & Ercolano (1200^^ 
which in turn are based on the methods presented bv iLucvl d 19991 ). 
These photon energy packets are collections of photons for which 
the total energy e remains constant, but the number of photons con¬ 
tained varies for different frequencies u. In the model, they are 
initiated at stars with frequencies selected randomly based on the 
emission spectrum of the star. The constant energy value e for each 
photon packet is simply the total energy emitted by star’s (luminos¬ 
ity L) during the duration At of the iteration divided by the total 
number of photon packets N: 


LAt 


(4) 


Photons that are emitted from a given ionizing source are prop¬ 
agated in random but isotropic directions. As soon as a photon 
packet is emitted, it will propagate for a path length I determined 
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by a randomly selected optical depth. The length I is determined by 
the position in which the next event will occur; either involving an 
interaction with the material after traversing a random optical depth 
given by 

r = -ln(l - r) ; r G [0,1), (5) 

(as detailed in iHarries & Howarthlll997h . or involving the crossing 
of a cell boundary. 

If the photon packet fails to escape a cell after traveling an 
optical depth r then its propagation ceases and an absorption event 
occurs: this is the situation when r > tcbll, tcell = k ■ dr 
where n is the opac ity term and dr the size of the cell (see also 
lErcolano et alj200^ . At this point, there are two possibilities based 
on the assumed effect of the diffuse field radiation; 


(i) Accounting for the diffuse field. The diffuse field refers to 
photons emitted by the recombination of electrons with ions. The 
effect of this radiation field is both to further ionize the gas (should 
the recombination photon be energetic enough) and to cool the gas 
(through extraction of energy by optically thin photons). The dif¬ 
fuse radiation field can be important for irradiating regions shad¬ 
owed by denser clumps along their line-of-sight with the ionizing 
source. Here, the total recombination coefficient (so-called ‘case- 
A’) is given by 


OLA 




( 6 ) 


where an{H°, T) is the recombination coefficient of hydr o gen a t 
temperature T in the quantum level n. IVemer & Ferlandl 1 19961) 
studied how a a changes by varying the temperature for the basic 
elements (H, He, Li, and Na). In TORUS-3DPDR, once an absorp¬ 
tion event occurs it it is followed by the immediate emission of a 
new photon packet from the same location and with the same en¬ 
ergy, but containing a different number of photons to account for 
different photon frequency. We also note that the diffuse compo¬ 
nent of radiation field can be created by the scattering of ionizing 
photons onto dust grains. Such tre atment has been performed e.g. 
lErcolano, Barlow. & Stored dlOOSh and we aim to include the dust 

contribution in a subsequent paper. _ 

(ii) Using the on-the-spot (OTS; lOsterbro^l 19891) approxima¬ 
tion (case-B) in which the diffuse field photons are assumed to con¬ 
tribute negligibly to the global ionization structure following ab¬ 
sorption. The recombination coefficient is then given by equation 


aB=aA-ai(IT°,r) = ^an(IT°,T) (7) 

n = 2 

This is justified in regions of simple geometry, for example where 
density gradients are small. In MC photoionization with the OTS 
approximation, once a photon packet is absorbed it is ignored and 
is assumed to either have been re-emitted with a frequency lower 
than that required for photoionization, or to provide a negligible 
further contribution to the ionization structure by causing further 
photoionization on only small scales. 

The energy density dU of a radiation field is given by 

dU = ^^do, (8) 

c 

where c is the speed of light and Ji, is the specific intensity and fre¬ 
quency v. A photon energy packet traversing a path I in a particular 
cell contributes an energy e{l/c)/At to the time-averaged energy 
density of that cell. Thus by summing over all paths I the energy 


density of a given cell (volume V) can be determined. Thus Eqn.[8] 
can be evaluated using the expression: 


Ttt Jix 
c 


dt' = 

cAt V 




(9) 


This is then used to obta in ionization fract ions by solving the ion¬ 
ization balance equation dOsterbrockll 1 9891) 

1 ^7tJ^a^{X')du 

niX^) a{Xi)nJ^^ hv ’ ^ ^ 

where n(X*), a(X*), a^(X’'), and are the number density 
of the ionization state of species X, recombination coefficient, 
absorption cross section, electron number density and the threshold 
frequency for ionization of species A* respectively. Using the MC 
estimators described in Eqn. Eqn.[To]can then be approximated 
as 


n(A'+i) _ e la^{X^) 

n(A») “ AtVa{Xi)n^ ^ hu 


This approach has the advantage that photon energy packets 
contribute to the estimate of the radiation field without having to 
undergo absorption events, thus even very optically thin regions 
are properly sampled. Photoionization calculations are performed 
iteratively, doubling the number of photon packets per iteration 
until the temperature and ionization fractions converge. Charge 
exchange recombination from O II, O III and N IV and ioniza- 
tion of N I and O I is in cluded following the prescription from 
iKingdon & Ferlan3 ( Il996h 


3.2 Chemistry in photoionized regions 


We include a range of atomic constituents: hydrogen, helium, car¬ 
bon, nitrogen, oxygen, neon and sulfur, for which we solve the ion¬ 
ization balance using Eqn. Note that the ionization balance is 
solved for every species in the calculation. The ionization states 
that we treat for metals are C (I-IV), N (I-IV), O (I-III), Ne (II-III), 
S (II-IV). The hydrogen, helium a nd C IV recombination rat es used 
by TORUS are calculated based on IVerner & Ferlan^ ) ll996l) . Other 
radiative recombination and dielec tronic recombination rates ar e 
calculated using fits to the result s of Nussbaumer & Storey! 1 19^ ), 
IPe gui^gnot, Petitiean, & BoissonI ( Il99lh or IShull & van Steenberd 
(Il982l) . The photoionization cross sections of all ato mic species in 
this pa per are calculated using the PHFIT2 routine from i Venter et al.l 
lll996h . 


3.3 Thermal Balance 


TORUS performs photoionization calculations that incorporate a 
range of atomic species and in which the thermal balance in each 
cell is calculated by iterating on the temperature until the heating 
and cooling rates match. Similarly to the photoionization calcula¬ 
tion (Eqn. [ID, the heating rate in a given cell is calculated based on 
the summation of trajectories of photon packets through the cell. 
This is used to estimate the heating c ontributions from the pho¬ 
toion ization of hydrogen and helium jWood. Mathis. & Ercoh^ 
2004) and the heating due to dust photoelectric ejection dLucv 
1999h . Adding up these three terms give the total heating rate. 

On the other hand, the cooling rate is initially calculated for 
the maximum and minimum allowed temperatures in the calcula¬ 
tion (3 X lO"* K and 10 K respectively by default in TORUS). This 
is then refined by bisection until the cooling rate matches the heat¬ 
ing rate. The cooling processes considered are free-free radiation. 
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hydrogen and helium recombination, dust cooling and collisional 
excitation of hydrogen and metals. 


4 COUPLING TORUS AND 3D-PDR 

The strategy we use to couple the two codes is to dismantle 3D-PDR 
and re-assemble it using the TORUS frameworlfl Although the it¬ 
eration procedure for determining the thermal balance, as well as 
the associated routines for estimating the heating and cooling func¬ 
tions, are identical to those of the 3D-PDR code, we have modi¬ 
fied the ray tracing scheme in order to overcome issues related to 
the large memory requirements. The most important feature how¬ 
ever, is the ability to use the MC radiation transport to compute 
an accurate structure of the UV field in any given arbitrarily com¬ 
plex distribution either for a single or multiple radiation sources. 
TORUS-3DPDR calculates steady states and its primary use is to 
post-process snapshots of hydrodynamical simulations. Below we 
discuss how the code calculates the UV field from the MC pho¬ 
toionization. In Appendix lAl we discuss the technical details fol¬ 
lowed to couple TORUS and 3D-PDR. 



Figure 1. This plot shows the extinction of the UV field in TORUS-3DPDR 
(blue crossed line) versus Eqn.[T3with the opacity term (green solid 
line) and without this term (red dashed line). TORUS-3DPDR reproduces 
precisely Ean. ll3l verifving that the UV estimate is accurate. 


4.1 Calculating the UV field from Monte Carlo 
photoionization 

In PDR modelling the UV field is usually prescribed at the ioniza¬ 
tion front and is then assumed to be attenuated according to some 
simple function (i.e. an exponential). In the coupled code we calcu¬ 
late the UV field using the propagation of photon packets in a pho¬ 
toionization calculation. The UV field in a given cell (in Draines) is 
the sum of the path lengths, I, of photon packets with frequencies 
in the UV (912 — 2400 A) and is given by 

where u is the mean radiation propagation v ector, Hp — 1.6 x 
10“® erg cm“^ is the Habing constant dHabind Il968h . and 
1.71 i s the conversi on factor between the Habing field to the Draine 
field dDrain ll978h . In practice we store the UV field as two vec¬ 
tors, one for positive and one for negative directions. This is impor¬ 
tant for multi-dimensional calculations, in particular if there are 
multiple sources (and the radiation field is not simply divergent 
about a point). 


5 CODE EVALUATION 

In this section we benchmark TORUS-3DPDR against 3D-PDR and 
MOCASSIN for simple one dimensional tests. The PDR abundances 
correspond to those of Milky Way solar undepleted abundances. 
The che mical network used he re is the most recent UMIST 2012 
netwo rk dMcElrov et al.l 1201 3l) in contrast to the I Woodall et al.l 
d2007h version used in B12. The cosmic-ray ionization rate is taken 
to be C(cr = 5 X 10“^^^ s“^. We also acc ount for the rate of H 2 for¬ 
matio n on grains using the treatment of ICazaux & Tieleli3 d2002L 
l2004h a nd we determine the dust tempera t ure fo llowing the treat¬ 
ment of lHollenbach. Takahashi. & TielensIdl99 ih in addition to the 
3D-PDR version presented in B12. 


^ Our future plans include a publicly available version of TORUS-3DPDR. 


5.1 Testing the UV intensity 


To check whether the UV estimator we use is accurate, we calcu¬ 
late the UV field as a function of distance from a test source. The 
source that we consider is the same as that used in the HII40 Lex¬ 
ington benchmark. We compare with our numerical estimate by in¬ 
tegrating over the UV band (912 — 2400 A) of the stellar spectrum 
and attenuate this value accounting for both the distance from the 
source and the gas opacity, i.e. 


X = 


Xo 


(13) 


where x is the attenuated Xo UV field, r is the distance from the 
source, and k is the opacity term (see Appendix O. We have com¬ 
pared the UV distribution (by invoking the on-the-spot approxima¬ 
tion) calculated by the TORUS-3DPDR with the above Eon. [Island 
we show the results in Eig. d] It can be seen that the agreement is 
excellent, verifying that our UV estimate is accurate. 


5.2 Consistency between 3d-pdr and torus-3dpdr 

In order to assess the accuracy of the PDR calculations after the 
coupling of the two codes, we perform a one-dimensional test to 
compare the coupled TORUS-3DPDR code against 3D-PDR. The 
setup of the test we adopt is v ery similar to the V2 model dis¬ 
cussed in the iRbllig et alj ( ITOOTL hereafter ‘R07’) paper. However, 
due to the recent updates of the 3D-PDR code described above, 
we do not aim at benchmarking TORUS-3DPDR against the various 
other one-dimensional PDR codes described in R07. The adopted 
model consists of a one-dimensional uniform density distribution 
with H-nucleus number density of nn = 10® cm“®. The size of 
the distribution is chosen so that the maximum visual extinction 
is Av.max = 10 mag. We assume that it is irradiated by a plane- 
parallel radiatio n field of stren gth xo = 10® multiples of the Draine 
radiation field dPrainell 19781) . and that the attenuation of the UV 
field is according to Eqn.[T] 

As described in B12, in order to emulate a semi-infinite slab 
we have considered cells aligned along two opposite HEALPIX rays 
at the i = 0 level of refinement while assuming very high optical 
depths in all other directions. Here, the 3D-PDR run uses A/av = 
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Table 1. Comparison between TORUS with the median value for all partic¬ 
ipating codes in the Lexington 2000 benchmark workshop for the standard 
H II region ( HII40) test. The values of the ‘median’ column are those listed 
in Table 4 of lErcolano et alj j2003h . The third column corresponds to the 
absolute emor 100 ■ |/toruS - ^med I/ 1^med I. where /med is the median 
value and the /toruS is the TORUS result for the intensity of each line. 


Line 

Median 

TORUS 

~EiTor% 

H/3/lO^^ ergs“^ 

2.05 

1.91 

7 

C III] 1907 + 1909 

0.070 

0.142 

51 

[N II] 122 ^lra 

0.034 

0.027 

26 

JN II] 6584 + 6548 

0.730 

0.801 

9 

JN II] 5755 

0.0054 

0.0077 

30 

iNlll] 57.3/im 

0.292 

0.263 

11 

JO l] 6300 -1- 6363 

0.0086 

0.0276 

69 

[0 II] 7320 + 7330 

0.029 

0.022 

32 

[0 II] 3726 + 3729 

2.03 

2.57 

21 

[0 III] 5007 -1- 4959 

2.18 

2.80 

22 

[0 III] 4363 

0.0037 

0.0108 

66 

[0 III] 52 + 88 fim 

2.28 

2.77 

18 

[Ne II] 12.8 fim 

0.195 

0.171 

14 

jNe III] 15.5 fim 

0.322 

0.350 

8 

jNe III] 3869 + 3968 

0.085 

0.110 

23 

IS II] 6716 -1- 6731 

0.147 

0.159 

8 

IS II] 4068 -1- 4076 

0.0080 

0.0062 

29 

IS III] 18.7 Atm 

0.577 

0.741 

22 

IS III] 9532 + 9069 

1.22 

1.23 

<1 


20 elements logarithmically distributed per Av dex with —5 ^ 
log(Av) ^ 1 in all cases, implying a total number of Afeiem = 
120 elements. The TORUS-3DPDR runs use a non-uniform mesh, 
clustered near the ionization front. The UV field is assumed to be 
plane-parallel impinging from one side while the attenuation of the 
field is calculated using Eqn.[T] 

Figure shows results for the benchmark model described 
above. We compare gas and dust temperature profiles, abundances 
of H I, H 2 , C II, C I and CO versus Ay, and local emissivities as 
well as emergent intensities (surface brightnesses) for some of the 
dominant cooling lines. The agreement between the two codes is 
excellent. We also note that we have performed all other tests de¬ 
scribed in the R07 work (i.e. VI, V3, and V4) and we have found 
excellent agreement as well. 


5.3 Photoionization Benchmarking 

The H II40 Lexington benchmark (lFerlandlll995l : Ipiguignot et al.l 
l200lh is a canonical test of photoionization codes. It involves mod¬ 
elling the ionization and temperature structure of a uniform den¬ 
sity medium around a massive star. It is usually performed in 
ID owing to the spheric ally symmetric nature of the problem. 
iHaworth & Harri^ ( 1201 2h demonstrated the veracity of the pho¬ 
toionization scheme in TORUS using this test, obtaining agreement 
with the ionization and temperature structu re computed by the well 
known CLOUDY code dFerland et alJl998ll . 

Table[T]compares TORUS with the median val ues of all partic¬ 
ipatin g codes in the HII40 Lexington benchmark dPeguignot et al.l 
l200lh . Overall we find good agreement between the lines emmisiv- 
ities. 


5.4 Photoionization - 1 - PDR test model 


In this section we perform a one-dimensional test application in 
which the integrated scheme is being tested to simulate the inter¬ 
action of a uniform-density cloud with ionizing radiation. The den¬ 
sity, riH, is taken to be 100 cm“® and the star is assumed to have a 
blackbody temperature of 4 x 10“^ K and radius 18.67R0 emitting 
photons from the left-hand edge of the grid. The grid is taken to 
be uniform and spherical consisting of 256 cells. These values cor¬ 
respond tothe HII40 simulation of the Lexington benchmarking 
dFerlandll995h . 

We convert the mean intensity into Draine units using the re¬ 
lation: 


X = 


(■2400 T ,, 
J912 

I.7IIT0 


(14) 


where J\ is the intensity of the radiation integrated over wave¬ 
length (here in A). Assuming monochromatic radiation at A = 
912 A, the integral of the above equation can be considered as a 
5-function and Eqn.[T4]can therefore take the form: 


_ A/uyC he 1 

^ 912A 1.71iTo ’ 


The above equation calculates the intensity of the UV radiation 
field (in Draine units) resulting from a star emitting A/uyC pho¬ 
tons per unit time placed at distance r. 

We perform two runs in which we explore the effect of the 
attenuation of radiation in the PDR scheme estimated i) directly 
from the MC approach as calculated in TORUS and described in 
1 14. H and ii) using the exponential relation given by Eqn.[T] TORUS- 
3dpdr selects as PDR cells any cell that has a temperature below 
3000K corresponding to an ionization fraction of Xi ^ 0.3. The up¬ 
per panel of Fig. shows the temperature profiles while the lower 
panel shows the fractional abundances of H II, H I, and H 2 . The 
ionizing radiation is impinging from the left hand side in the above 
diagrams (the position of the star defines the centre of the Carte¬ 
sian co-ordinate system). In addition, in the upper panel we o ver- 
plot the results obtained from MOCASSIN dFrcolano et al.l200^ for 
the calculations concerning the photoionized region, while we also 
overplot those obtained from 3D-PDR for the photodissociation re¬ 
gion. We find that the integrated code TORUS-3DPDR is in very 
good agreement with MOCASSIN. For the particular case of the 
MC-attenuation, we have performed an additional 3D-PDR calcu¬ 
lation in which we have included the additional extinction due to 
the distance of the point source from each PDR cell as described 
by Fqn.[^ As is shown in the upper panel of Fig.[^ the PDR tem¬ 
perature profile obtained from the MC TORUS-3DPDR technique is 
reproducible only when the distance from the source is taken into 
account. Note, however, that H 2 and H I are not strongly sensitive 
to that. 


6 APPLICATIONS 

In this section we present the capabilities of the code in simulating 
three-dimensional density structures interacting with one- or two- 
sources, either by invoking the on-the-spot approximation or in¬ 
cluding the diffusive component of the ionizing radiation, thus ex¬ 
ploring the importance of these components in regulating the chem¬ 
istry of PDRs. To do this we consider an oblate spheroidal cloud 
with a uniform density distribution. The choice of such an object 
was based on the assumption it is not spherically symmetric hence 
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Figure 2. Comparison between 3D-PDR (solid lines) and TORUS-3DPDR (filled circles) for the benchmark model V2. Panel (a) shows the gas and dust 
temperature profiles. Panel (b) shows the H and H 2 abundances while panel (c) shows the C II, C I and CO abundances versus the visual extinction Ay. Panel 
(d) shows the local emissivities for [C I] 370/rm and [O I] 63^m and panel (e) shows the local emissivities for three different CO line transitions. Finally, panel 
(f) shows the surface brightnesses for the dominant cooling lines. In all cases the agreement is excellent. The discrepancy observed in the bottom left panel at 
Ay ~ 5 is due to noise resulting from the grid. 


not reproducible by an equivalent one-dimensional approach, while 
at the same time it offers simplified results. 

We consider a cubic computational domain with size R — 
12 pc consisting of 64® cells of equal volume. Figure |4] shows the 
simulation setup. The center of the spheroid defines the center of 
a Cartesian co-ordinate system. The semi-major axis is taken to 
be i?i = 4 pc, the semi-minor axis is taken to be R 2 = 2 pc 
and the density is assumed to be tie = 1000 cm“®. This means 


that the visual extinction, ^y, can reach values up to ~ 8 mag 
along the y—axis. We can consider this structure as a cloud of mass 
Me = 1655 Mq. The spheroid is embedded in an ambient medium 
of uniform density tia = lcm~®. In all simulations we use the 
MC photoionization routine described in 113. H and we use an£ = 0 
HEALPIX level of refinement as explained in 116.11 

We perform three simulations. In the first simulation (SI), we 
place a single exciting source at a distance of D = 4.86 pc from 
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Figure 3. Upper panel: Gas temperatures obtained by TORUS-3DPDR 
(black solid line) and 3D-PDR (red dashed, and blue dot-dashed). The blue 
dot-dashed line coiTesponds to the attenuation of the UV field as described 
by Eqn. ^ whereas the red dashed line is as described by Eqn. [3] Black 
crosses (-H) show a comparison with the MOCASSIN code. Lower panel: 
fractional abundances of hydrogen for different phases. The black solid line 
is ionized H II, the green solid line is atomic H I, and the blue-gray solid 
line is molecular H 2 . The green dot-dashed line is atomic hydrogen ob¬ 
tained without PDR calculations (i.e. as in the Lexington test). All these 
correspond to the TORUS-3DPDR code. The green and blue dashed lines 
correspond to 3D-PDR runs using Eqn.^to account for the attenuation of 
the radiation field. 


the center of the spheroid emitting A/LyC = 1.2 x 10^*^ photons per 
second and invoking the on-the-spot approximation. In the second 
simulation (S2), we use the same setup but we switch on the con¬ 
tribution of the diffuse component of the radiation field. Thus by 
comparing simulations SI and S2 we examine the effect of the dif¬ 
fuse radiation field on the abundance distribution. In the third sim¬ 
ulation (S3) we replace the single exciting source with two sources 
emitting half the number of photons per second, compared to the 
SI and S2 cases. In particular, we place the two stars at positions 
(x, y) = (±1.98, 4.44) pc which ensures that the distance D from 
the center of the spheroid is kept the same as in S1 and S2. Each of 
the two stars emit Ni^yc = 6 x 10'^® photons per second and the 
diffuse component of radiation is taken into account. Therefore, in 
the S3 simulation we test the capabilities of the code in treating 


S] 

S3i © * 

\ 

,S2 

* pS3, 

/ 

/ 

D /' 



r2 

= 1 cm~^ 

- 


-6 -4-20 2 4 6 

X (pc) 


Figure 4. Setup of the thi'ee-dimensional oblate spheroidal cloud case 
discussed in ^ The x— and y— axis dimensions are in pc. The shaded 
spheroid has dimensions of Hi =4 pc and R 2 = 2 pc and it has a uniform 
density with ng = 10^ cm“^. The center of the cloud defines the center 
of a Cartesian co-ordinate system. The spheroid is embedded in an ambient 
medium with density np^ = 1 cm“^. The middle black point corresponds 
to the position of the ionizing source for the simulations, SI m6.2\ and S2 
(gO). The gray points, S3i and S32, correspond to the two positions of the 
stars in the S3 simulation (gM). All stars are placed at a distance D from 
the center of the oblate spheroidal cloud. 


PDRs interacting with multiple sources and we examine the conse¬ 
quent effects introduced in the abundances distribution as well as 
on the distribution of the local emissivities of the most important 
cooling lines. 

The stellar spectra in the models presented here are blackbod- 
ies, which are probabilistically randomly sampled by the MC pho¬ 
ton packets. There are therefore occasionally high frequency ion¬ 
izing photon packets, that can propagate far into the neutral gas 
since the photoionization cross section is proportional to the in¬ 
verse cube of the photon frequency, i.e. ex (zztrsh/t^)^ where the 
zztrsh is the ionization threshold (e.g. 13.6 eV/h for hydrogen). 
These high frequency packets are responsible for regions of local¬ 
ized heating exterior to the main ionization front. In the case of our 
three-dimensional models, the calculations with a single star con¬ 
sider a source that has 7? = 18 Rq radius and Ten = 5 x 10"* K 
effective temperature (whereas in the multiple sources simulation 
each source has R = 12.7 Rq with half that effective tempera¬ 
ture), and so have a higher probability of emitting high frequency 
photon packets. This is why the temperature distribution within the 
ellipse and particularly in the ionization front is subject to spurious 
heating which has a significant (but very local) effect in determin¬ 
ing the abundances distribution. These areas can be seen in almost 
all panels of Fig.|7]and Fig.l^and are located close to the a; = 0 pc 
axis (i.e. x ^ —1 pc and a; ~ 1 pc) and with y > 1 pc. 

This can be alleviated by using more photon packets (so the 
energy per packet is reduced) however at the cost of computational 
expense. We will thus exclude this effect from all discussion be¬ 
low. We also note that in the particular case of simulation S3 and 
because the ionizing photons have lower frequency than in the S1 
and S2 cases, this effect is not present. 
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Figure 5. Cross section plots of the SI simulation using different levels of HEALPIX refinement. Top row is gas temperature and bottom row is the abundance 
distribution of H 2 . From left to right £ = 0 (12-rays), £ = 1 (48-rays), 1 = 2 (192-rays) of refinement. It can be seen that there are no significant differences 
for the tests presented here when moving to high angular resolution. For all simulations presented in this paper we shall use 1 = 0. The gray color in the 
background of the bottom row indicates a zero value. 


6.1 Angular resolution in the PDR regime 

As a first test we explore the effect of using different levels of 
HEALPIX refinement corresponding to different angular resolutions 
per direction in each cell of the computational domain. We perform 
two additional simulations to SI, one with £ = 0 (correspond¬ 
ing to 12 HEALPIX rays) and one with £ = 2 (corresponding to 
192 HEALPIX rays). Figure shows the changes observed for the 
gas temperature and H 2 abundance distribution. Overall we find 
excellent agreement in the gas temperature distribution and minor 
changes in the H 2 abundance distribution, showing that even at low 
angular resolution (i.e. at f = 0) our results remain consistent with 
those at higher angular resolution (i.e. I = 2). Given also that the 
computational expense increases approximately 4 times by increas¬ 
ing the HEALPIX level (since there are 4 times more rays emanated 
from each cell; see Appendix |A3] for further discussion) in this pa¬ 


per we will use £ = 0 for all simulations. We note that in more 
complicated simulations with non-uniform density distribution, a 
higher level of angular resolution must be used. 


6.2 Single source and on-the-spot approximation (SI) 

Figure|6]shows the outcomes of the S1 simulation. We show cross- 
section plots (at z = 0 pc plane) of the gas temperature (Fig.|6]4), 
the abundance distributions of FI 2 , FI I, C H and C I (Fig.|^-E re¬ 
spectively) and the local emissivity of [C I] 609 fim (Fig.[^). The 
ionizing source interacts with the rarefied medium resulting in an 
H H region. The dense oblate spheroidal cloud on the other hand 
shields the propagation of the UV photons behind it. The resultant 
temperature in the FI H region is ~ 10^ K while in the shadowed re¬ 
gion it is ~ 50 K, and in the innermost part of the spheroid ~ 10 K. 
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This extra heating, that shifts the temperature from 10 K to 40 K 
in the ambient medium behind the spheroid, is due to the interaction 
of that rarefied gas with cosmic rays, and which do not attenuate in¬ 
side this computational domain. Note that the gas temperature and 
ionization fraction transition from the ionized medium to the PDR 
is very abrupt. 

The shielding of the UV radiation allows the formation of 
atoms and molecules in the interior of the cloud. In the modelled 
cloud, this shielding extinguishes the UV radiation to such an ex¬ 
tent that is able to photodissociate H 2 to form H I while it cannot 
photoionize atomic hydrogen. This creates a layer rich in atomic 
hydrogen as shown in Fig.|6fl. Here, the north surface shows an 
increase in the abundance of atomic hydrogen. Immediately behind 
the ionization front (i.e. towards the interior of the cloud), the free 
hydrogen atoms transition to molecular hydrogen as shown in Fig. 
[6jl. Figuresl^ and E show the distribution of C II and C I respec¬ 
tively. By comparing these two panels it can be seen how the C II to 
C I transition occurs in the spheroid. Figure|^ plots the local emis- 
sivity of [C I] 609 ^m. This line is predominantly emitted from a 
thin region located at low optical depths i.e. close to the ionization 
front, and its peak is at Av 1 mag. 

We emphasize that the H I-to-H 2 transition is known to oc¬ 
cur in a very thin zone of the PDR and consequently in order to 
properly resolve it we need to increase the resolution of the grid. In 
the work presented here, none of the 3D simulations use high spa- 
ti al resolution in orde r to study this transition in detail. As shown 
in lOffner et alj (l2013h . such a low spatial resolution (i.e. not being 
able to resolve the H I to H 2 transition) does not change the overall 
chemistry of the PDR. However, by using the AMR techniques as 
described in Appendix lAl TORUS-3DPDR is able to perform such 
detailed studies at the cost of computational expense. 


6.3 Effect of the diffuse radiation (S2) 


To explore the effect of diffuse radiation on the abundance distribu¬ 
tion of different species, we repeat the SI simulation in which we 
are additionally considering a full treatment of the UV field. The 
diffuse component of the radiation field penetrates further into the 
cloud, thus increasing its temperature and changing its chemistry 
and line emission. We focus here on the changes observed in the 
oblate spheroidal cloud itself between the SI and S2 simulations, 
hence we will not include in our discussion the ambient medium. 
To show the effects of the interaction of the diffuse field we plot in 
the left column of Fig. Across sections of the spheroid, coloured 
according to the percentage of the relative change, ct%, between 
the two simulations, given by the relation 


a% 


-Sj 


Sj 


( 16 ) 


where Si and Sj are the values of the i and j simulation that we 
compare. 

The relative change in the gas temperature is shown in Fig. 
EK. Here, although the temperature remains the same close to the 
ionization front (with a relative change of < 20%), it gets hotter in 
the inner part of the cloud (with a relative change of ~ 90%). In 
particular, at high optical depths we find that the spheroid in the S1 
simulation has an average temperature of Tg^s ~ 10 K, whereas in 
the S2 simulation it has Tgas ~ 100 — 200 K. As a consequence, the 
abundance of atomic hydrogen is higher in the S2 simulation reduc¬ 
ing at the same time the abundance of H 2 . The H I-to-H 2 transition 
has now been shifted further inside the cloud (but still at low optical 
depths; see Fig.[7j3). Similarly, C I is primarily emitted by a region 


closer to the centre of the ellipsoid. Since the peak of C I occurs 
in general in a quite thin zone of the PDR, such a shift may cause 
significant changes in its abundance distribution through the ellip¬ 
soid as observed in Fig.EJi. However, in the innermost region, both 
SI and S2 simulations are in excellent agreement. Hence inclusion 
of the diffuse component of the radiation field primarily affects the 
abundance distribution at low optical depths, while it is able to heat 
up the gas at intermediate optical depths. 

Figure [8] shows cross section plots of the ellipsoid where the 
local emissivity of CO J = 1 — 0 and CO J = 7 — 6 transi¬ 
tion lines are mapped (left and middle panels for the SI and S2 
simulations respectively). The difference in gas temperature at the 
inner part of the cloud due to the diffuse component of radiation 
leads to discrepancies in the emissivities of both J transition lines 
with the high-J to be up to ~ 5 orders of magnitude brighter in 
the S2 simulation. This shows that high-J CO lines may be used 
as gas temperature diagnostics in regions of intermediate optical 
depths where the diffuse radiation can penetrate the cloud deeper. 
It also shows that a realistic treatment of the UV field is neces¬ 
sary when comparing observational with numerical data such as i.e . 


of the Orion Bar (i.e.lTielens et al.||l993 

; Ivan der Werf et al.ll 19961; 

lAndree-Labsch. Ossenkonf. RbllidiTOM 

). 


6.4 Effect of multiple sources (S3) 

The effect of multiple sources interacting with PDRs is particularly 
interesting and very poorly studied due to the lack of numerical 
codes that are able to treat such systems. Here, we discuss the ca¬ 
pabilities of TORUS-3DPDR in modelling such PDRs by revisiting 
the S2 simulation described above, in which we replace the ioniz¬ 
ing source of A/LyC = 1-2 x 10®° with two identical sources 
of A/Lyc = 6 X lO"*® s“^. Both sources in this simulation have the 
same distance D from the center of the oblate spheroidal cloud as 
in the SI and S2 simulations. The diffuse component of the radia¬ 
tion field is taken into consideration, hence this setup ensures that 
any changes observed between S2 and S3 will be a result of the 
structure of the radiation field emitted by two sources instead of a 
single source. 

The right panel of Fig. |7] shows the relative change between 
S2 and S3 simulations for the gas temperature, the atomic hydro¬ 
gen abundance, and the atomic carbon abundance. We find that the 
gas temperature shows fluctuations of the order of < 30% for the 
entire cloud, corresponding to ~ 20 K (see Fig. E^). The spatial 
distribution of the UV field has an impact in the abundances distri¬ 
bution between S2 and S3 particularly at low optical depths. This 
can be seen in Fig.EP andEP, where we compare H I and C I. Simi¬ 
larly to what has been discussed in il6.3l the UV distribution is now 
less strong as in S2 resulting to a shift of the H I-to-H 2 transition to 
lower optical depths, thus also increasing H 2 . Such a behaviour is 
also observed in C I. 

The middle and right columns of Fig. compares CO J = 
1 —OandCO J = 7 —6 for the S2 and S3 simulations. In both cases 
the differences observed are not significant in comparison with S1 
and S2, however at the inner part of the cloud the J = 7 — 6 line 
emissivity is higher by a factor of ~ 2 reflecting the gas temper¬ 
ature difference of ~ 20 K between S2 and S3. We argue that in 
more complicated structures involving multiple sources the effects 
introduced by the realistic UV treatment can be of significant im¬ 
portance. Furthermore, we plot in Fig.|^the relative changes in the 
line emissivities for the S2 and S3 simulations of C II 158 /rm, C I 
609 lira, O I 146 lira and CO J = 1 — 0. In general, C II 158 
remains in principle unaffected when the UV radiation field struc- 















T0RUS-3DPDR 11 




x[pc] 


x[pc] 


x[pc] 


10 -' 10 -“ 10 -“ 10 -* 
Fractional abundance 


10 -* 10 -“ 10 -“ 10 -* 
Fractional abundance 


10-27 


10-2S 

Ij. [erg cm-“ s-‘ ] 


10 - 


Figure 6. Cross section plots for the SI simulation. The axes are in pc. Panel A shows the gas temperature distribution. Panels B-E show the distribution of 
the fractional abundances of molecular hydrogen (H 2 ), atomic hydrogen (H I), ionized carbon (C II), and atomic carbon (C I) respectively. The abundances 
of these species ai'e with respect to the total hydrogen nucleus abundance. Panel F shows the local emissivity of [C I] 609/rm. The colour bar in panel F is in 
units of erg cm”® s”In panels B-F the gray background indicates a zero value. 


ture is taken into account, however CI 609 nm and in particular CO 
J = 1 — 0 show strong dependency under a multiple source treat¬ 
ment. These point to the conclusion that emission lines do become 
dependent on the spatial structure of the radiation field even when 
the total intensity received is the same as the intensity emitted by 
a single source at the same distance. This is in agree ment with the 
results presented in B12 as well as with the findings of lOffner et al.l 
ll2013h . Thus three-dimensional codes such as TORUS-3DPDR are 
often needed in order to compare as closely as possible numerical 
simulations with observations. 


7 CONCLUSIONS 

In this paper we present TORUS-3DPDR, a three-dimensional astro- 
chemistry code which is able to treat simultaneously photoioniza¬ 
tion and photodissociation regions with arbitrary geometrical and 


density distributions interacting with single or multiple sources. 
The code uses a Monte Carlo-based scheme to calculate the prop¬ 
agation of the ionizing radiation including its diffusive component 
and a HEALPIX-based ray-tracing scheme to calculate the column 
densities in each cell and along any direction on the celestial sphere. 
We have performed one-dimensional and three-dimensional simu¬ 
lations demonstrating the capabilities of the code. 

We performed a one-dimensional PDR calculati on to explore 
the c onsistency of TORUS-3DPDR against 3D-PDR ( iBisbas et al.l 
l2012h and we found excellent agreement between the two codes 
in reproducing the temperature profile as well as the abundances of 
species and the local emissivities of the major cooling lines. The 
simulation setu p was identical to a standard PDR benchmarking 
test discussed in iRollig et alj||2007h . although our codes include the 
updated UMIST network as well as an updated H 2 formation heat¬ 
ing mechanism and thus we can not compar e TORUS-3DPDR di - 
rectly against the numerous codes presented in iRollig et al.l ( l2007t) . 











































12 


T. G. Bisbas et al. 




10 20 30 40 50 60 70 80 90 100 

CT% 



cr% 


Figure 7. Cross section plots showing the relative change between SI and S2 simulations (left column) and S2 and S3 simulations (right column). Thus the left 
column corresponds to the effects observed when including the diffuse component of the radiation field while the right column corresponds to when a multiple 
source treatment is taken into consideration. From top to bottom, we plot the relative change of the gas temperature, atomic hydrogen abundance, and atomic 
carbon abundance. The noise observed at the innermost part of the spheroid in Panel B (and which corresponds to a 20 K temperature difference) is due to 
the thermal balance accuracy used. The x— and y— axes are in pc. Relative changes < 5% and > 100% are not shown. As a result of the additional diffuse 
component of the radiation field in S2, the ellipsoid is warmer (by 40 K on average) thus changing the abundances of species accordingly and paiticularly 
at low optical depths (left column). On the other hand, a multiple source treatment (right column) shows that although the gas temperature remains remarkably 
similar between the S2 and S3 simulations, H I and C I show a dependence on the UV field structure at low optical depths. 
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Figure 8. Cross section plots showing local emissivities of CO for the transitions J = 1 — 0 (upper row) and J = 7 — 6 (lower row). The first column 
corresponds to the SI simulation, the second column corresponds to the S2 simulation and the third one to the S3. The units of the logarithmic colour bar 
are in erg cm“® s“^. The diffuse component of the radiation field heats up the inner part of the cloud by a few K which result in changing the emissivity of 
the different J transition lines with the high ones to be more affected, as we observe by comparing the SI and S2 simulations. However, a multiple source 
treatment (simulation S3) results in an emissivity map not significantly different from that of the S2 simulation, although it affects the emission at low optical 
depths (i.e. close to the ionization front). 


We then advanced our benchmarking calculation by per¬ 
forming a one-dimensional calculation including both the MC 
photoionization treatment as well as the PDR treatment in a 
uniform density cloud. The photoionization calculations repro¬ 
ducing the H II region followed the Lexington HII40 bench¬ 
marking case and showed excellen t agreement with the MO¬ 
CASS I N code ter colano et"^ I 2 OO 3 I: IWood, Mathis. & Ercolanol 


|2004 lErcolano, Barlow. & Storey 2 OO 5 I) . Once the value of the 
ionization fraction of hydrogen fell below Xi < 0.3, correspond¬ 
ing to a gas temperature of T 3000 K, the PDR calculations 
were switched on. The thermal balance obtained in the PDR using 
TORUS-3DPDR was in agreement with 3D-PDR only when we in¬ 
cluded the extinction due to the distance of each PDR cell for the 
ionizing radiation emitted by the ionizing point source. We thus ar¬ 
gue that this factor is important in order to obtain more accurate 
abundance distributions in the PDR and particularly a more accu¬ 
rate gas temperature profile. 

We have demonstrated the three-dimensional capabilities of 
TORUS-3DPDR by treating an application in which we considered 
an oblate spheroidal cloud irradiated externally by either a sin¬ 
gle source or two sources. Although this simulation setup is sim¬ 
plified in three-dimensional space, it is not reproducible by one¬ 
dimensional codes. We have performed control runs in which we 
have invoked the on-the-spot approximation and additionally the 
effects introduced by the diffuse component of radiation field. Our 
results for this application can be summarized as follows: 


(i) The effects introduced by the diffuse component of the radia¬ 
tion field in the single ionizing source case (simulations SI and S2), 
have been found to affect both the outer parts (low optical depths) 
and the inner parts (high optical depth) of the cloud. We found that 
the gas temperature has been increased as a result of the diffuse ra¬ 


diation by an average of ~ 40 K and by ~ 100—200 K individually 
at high optical depths. This in turn resulted in a higher abundance of 
H I while also increasing locally the abundance of C II compared to 
what we found when invoking the on-the-spot approximation (SI). 

(ii) The effects introduced by the interaction of two sources in¬ 
stead of a single source (simulations S2 and S3) can heat up the 
gas at the inner part of the cloud by 20 K which results to a 
~ 2 times higher line emissivity of the CO J = 7 — 6 transition. 
We argue that the impact of multiple source treatment is potentially 
higher at even higher CO J transitions. In addition, the abundance 
distributions of the species examined are strongly correlated with 
the geometric distribution of the radiation field. 

(iii) The local emissivities of the fine-structure lines which are 
predominantly emitted at low optical depths, depend on the spa¬ 
tial shape of the radiation field (simulations S2 and S3). The local 
emissivities of CO transitions are strongly affected, showing that 
we require realistic and detailed UV structures when comparing 
models with observational data. 

In addition, we have performed a test in which we increased 
the angular resolution (HEALPIX level of refinement) in the SI sim¬ 
ulation and examined the gas temperature structure and H 2 abun¬ 
dance distribution. We found that even at the lowest possible reso¬ 
lution (£ = 0 which also compromises a low computational cost) 
TORUS-3DPDR provides acceptable results for the simulations pre¬ 
sented in this paper. For more complicated density distributions 
however, a higher level of HEALPIX refinement is needed though 
at the cost of computational expense. 

While most existing PDR codes consider one-dimensional 
density profiles and very simplified radiation fields, TORUS-3DPDR 
brings us a step closer to simulating realistic three-dimensional 
H II/PDR complexes interacting with multiple sources. Future ver- 
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Figure 9. Cross section plots showing the relative change of the local emissivities for various cooling lines for the S2 and S3 simulations and where in both 
cases diffuse radiation has been taken into consideration. Both axes are in pc. Panel A shows the relative change for the [C II] 158/rm line. Panel B for the 
[C I] 609/rm, panel C for [O I] 146/rm and panel D for CO (1-0). From these plots we see that overall multiple sources have an impact on the emissivities 
of the cloud, although in both S2 and S3 simulations the total intensity of UV radiation used was the same. As explained also in Fig|7] the areas close to 
ionization front and to the a; = 0 pc axis of symmetry are subject to spurious heating due to a small number of high energy photons. 


sions of TORUS-3DPDR will include treatments of dust radiation 
transfer and X-ray chemistry in the photoionization and photodis¬ 
sociation regions. 
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APPENDIX A: TECHNICAL CHARACTERISTICS 

AI Ray tracing on a domain decomposed grid 

TORUS-3DPDR uses the octree Adaptive Mesh Refinement (AMR) 
grid native to TORUS, including domain decomposition. Domain 
decomposition splits the grid into subsets, each of which are han¬ 
dled by processors with independent private memory. This provides 
the advantage of alleviating intensive memory requirement of of 
3D-PDR. However, during ray tracing, rays must traverse cells on 
the grid that belong to different domains to that from which the ray 
originate, in which the thread doing the ray tracing has no knowl¬ 
edge of the conditions. To overcome this issue. Message Passing 
Interface (MPI) communication is required between the thread do¬ 
ing the ray tracing and the thread through which the ray is passing. 
To facilitate this, we iterate over domains, with one performing ray 
tracing operations and the others acting as servers. MPI does not 
therefore offer significant speedup in the ray tracing component 
of the calculations. It does however offer speed up in other com¬ 
ponents of the PDR calculation and offers excellent spe edup for 
the M onte Carlo radiative transfer part of the calculation jHarriesI 

[20TI) . 


approximately 22 CPU minutes. The PDR calculations in the three- 
dimensional (excluding the Monte Carlo photoionization part) re¬ 
quired 148 CPU hours for ^ = 0 HEALPIX refinement, ~ 522 
CPU hours for ^ = 1 and ~ 1251 CPU hours for 1 = 2. These 
show that an £ = 2 simulation is about an order of magnitude 
more expensive than one at ^ = 0. We are currently upgra ding the 
OpenM P parallelization and with recent developments bv lHarriea 
( I 2 OI 5 I) . the computational cost will be reduced in future versions 
of TORUS-3DPDR. These further optimizations will make compu¬ 
tationally more tractable complicated post processing of 3D simu¬ 
lations. 


APPENDIX B: FLOWCHART OF TORUS-3DPDR 

Figure im shows the basic flowchart describing how TORUS- 
3dpdr converges in evaluating the photoionization (upper half of 
the Figure) and PDR (lower half of the Figure) calculations. Each 
solid box represents a DO-loop over grid cells. The dashed box rep¬ 
resents iterations over the ionization and the thermal balance in the 
photoionization calculations. 

The code starts by reading all input model parameters, such 
as the density distribution, the abundances of species, and the loca¬ 
tion and properties of each ionizing source. It then starts the Monte 
Carlo based propagation of ionizing photons and it iterates in or¬ 
der to obtain equilibrium in the H It region. Once TORUS-3DPDR 
reaches photoionization convergence, it stores the direction of the 
UV radiation field in each cell. The code then proceeds further 
in performing the PDR calculations. In particular HEALPIX-based 
rays emanate from each cell flagged as PDR, and it moves along 
them to calculate the column densities for every cell in every direc¬ 
tion. A cell is flagged as PDR if the value of the hydrogen ionization 
fraction is below Xi < 0.3, corresponding to a gas temperature of 
T ~ 3000 K. 

Any further PDR calculation fo llows the flowchart presented 
in Appendix A of iBisbas et alJi2012l) . Once TORUS-3DPDRis fully 
converged it writes the outputs and terminates. 


A2 On-the-fly calculation 

As explained in iBisbas et alj ( l2012h . the 3D-PDR code stores in¬ 
formation on every intersection point for every ray. Although this 
offers the ability to perform PDR calculations directly on any SPH 
or grid-based simulation without further modifications, it has dis¬ 
advantages with regards to the large memory requirements. In the 
new coupled code, we perform calculations on a cell-by-cell ba¬ 
sis. We take advantage of the fact that TORUS moves recursively 
through the AMR grid throughout the calculations. 


A3 Hybrid parallelization 

In addition to the domain decomposition we simultaneously par¬ 
allelize the code using Open Multi-Processing (OpenMP) over in¬ 
tensive loops in the calculation. On each computational node, one 
(or more) cores are assigned the role of MPI threads and the oth¬ 
ers are used for OpenMP parallelization. Use of MPI and OpenMP 
in this hybrid manner provides us with a flexible parallelization 
scheme that both reduces calculation wall time and memory de¬ 
mands. For the tests performed here, we find that the computa¬ 
tional expense depends strongly on each different application. The 
one-dimensional HII/PDR calculation presented in 115.41 required 


APPENDIX C: GEOMETRICAL DILUTION FACTOR 


The flux, F, at a given point point R is given by the equation 


F = 


Lq 


(Cl) 


where Lq is the luminosity of the source (in this case a single star). 
At the position of the ionization front (IF), the above equation leads 
to 


Lq = 4:71 Rip Fif. (C2) 

At a given position, r, inside the PDR (hence at a distance R = 
Rip + r from the star), Eqn. lEB takes the form 


Fr = 


Lq 


47r(i?iF + rY ’ 
which when combined with Eqn.([C2j leads to 


F) = Fif 


Rl 


(Fif -f r)2 ■ 


(C3) 


(C4) 


From the above it turns that the factor Fif/(Fif + corre¬ 
sponds to the geometrical dilution of the UV radiation when ac¬ 
count is taken of the dilution inside the H II region. Hence, if we 
additionally take into account the exctinction of the radiation along 
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Figure Bl. Flowchart of TORUS-3DPDR 


distance r, it leads to Eqn. 0). Neglecting the dilution of the ra¬ 
diation inside the H II region, the corresponding factor is simply 
(47rr^)“^, leading to Egn. jlSb . We note that a full treatment of 
the H II/PDR complex using the MC method takes into account 
the dilution inside the ionized region, hence it is in agreement 
with Eqn.lO. In addition, for simplified cases such as in the plane- 
parallel approx imation, TORUS-3D PDR is able to treat the UV field 
as discussed in iRollig et all ( l2007h . 































